##########################################################################################

library(data.table)
library(optparse)
library(ArchR)
library(dplyr)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_all_file"), type = "character"),
    make_option(c("--type"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
  ## 所有的细胞的,计算maxdelt
  comine_data_all_file <- "~/20231121_singleMuti/results/qc_atac_v2/germ/testis_combined_peak.combineRNA.qc.Rdata"

  ## 输出
  out_path <- paste0("~/20231121_singleMuti/results/tf_regulators_", cor , "_" , maxDelta , "/germ/")

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_all_file <- opt$comine_data_all_file
out_path <- opt$out_path
type <- opt$type

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据

b <- load(comine_data_all_file)
## testis_combined_peak_combineRNA


##########################################################################################
# Identify Correlated TF Motifs and TF Gene Score/Expression
##########################################################################################
# To identify 'Positive TF regulators', i.e. TFs whose motif accessibility 
# is correlated with with their own gene activity (either by gene score or gene expression)
seGroupMotif <- getGroupSE(ArchRProj=testis_combined_peak_combineRNA, useMatrix="MotifMatrix",groupBy="cell_type")
seZ <- seGroupMotif[rowData(seGroupMotif)$seqnames=="z",] # Subset to just deviation z-scores

# identify the delta in z-score between all cells
res <- c()
for( i in 1:nrow(assay(seZ))){
  motif_name <- seGroupMotif@elementMetadata$name[i]
  cell_name <- names(which.max(assay(seZ)[i,]))
  zscore <- max(assay(seZ)[i,])
  delta <- zscore - min(assay(seZ)[i,])
  tmp <- data.frame( motif = motif_name , cell_type = cell_name , zscore = zscore , delta = delta )
  res <- rbind( res , tmp )
}

out_file <- paste0( out_path , "/Motif_maxScore." , type , ".tsv")
write.table( res , out_file , row.names = F , sep = "\t" , quote = F )

## 原始的开放程度列表
score_matrix <- data.frame(assay(seZ))
score_matrix$motif <- seGroupMotif@elementMetadata$name[1:nrow(score_matrix)]
out_file <- paste0( out_path , "/Motif_Score." , type , ".tsv")
write.table( score_matrix , out_file , row.names = F , sep = "\t" , quote = F )
